回帰 (Regression)
以下のYouTube動画を視聴して、回帰 (Regression)の概念を理解してください。動画は英語ですが、日本語字幕機能をオンにすることもできます。
Please watch the following youtube video to get the idea of regression
analysis
実証例のデ一タ可視化 (Data Visualization)
研究課題:少人数教育が学生の成績に与える定量的効果は何か?
以下のデータセットCASchoolsを使用します。
ここで使用するデータは、カリフォルニア州の420のK-6とK-8の全地区のデータで、1998年と1999年のデータが利用可能です。テストのスコアは、5年生を対象に実施されたスタンフォード9標準化テストのものです。
学校の特徴(地区全体の平均値)には、入学者数、教師数(「フルタイム換算」で測定)、教室あたりのコンピュータの数、生徒一人あたりの支出が含まれます。
生徒の人口統計学的変数は、地区全体で平均化されています。
人口に関連する変数には、生活保護プログラムであるCalWorks(旧AFDC)の生徒の割合、割引価格の昼食の対象となる生徒の割合、英語学習者(英語が第二言語である生徒)の割合が含まれています。
We are going to use the following dataset CASchools. The
data used here are from all 420 K-6 and K-8 districts in California with
data available for 1998 and 1999. Test scores are on
the Stanford 9 standardized test administered to 5th grade students.
School characteristics (averaged across the district)
include enrollment, number of teachers (measured as “full-time
equivalents”, number of computers per classroom, and expenditures per
student. Demographic variables for the students are
averaged across the district. The demographic variables include the
percentage of students in the public assistance program CalWorks
(formerly AFDC), the percentage of students that qualify for a reduced
price lunch, and the percentage of students that are English learners
(that is, students for whom English is a second language).
# データはAERというパッケージに含まれています
# 下記のパッケ一ジをlibrary()関数で読み出します
library(AER)
# 利用したいパッケ一ジがなかったら、
# 下記のようにパッケ一ジをインスト一ルしておきます (初回のみ); 例として:
# install.packages("AER")
# デ一タを読み込む
# AER::CASchools
# load the data set onto the workspace
data("CASchools") # dataframe structure
str(CASchools)'data.frame': 420 obs. of 14 variables:
$ district : chr "75119" "61499" "61549" "61457" ...
$ school : chr "Sunol Glen Unified" "Manzanita Elementary" "Thermalito Union Elementary" "Golden Feather Union Elementary" ...
$ county : Factor w/ 45 levels "Alameda","Butte",..: 1 2 2 2 2 6 29 11 6 25 ...
$ grades : Factor w/ 2 levels "KK-06","KK-08": 2 2 2 2 2 2 2 2 2 1 ...
$ students : num 195 240 1550 243 1335 ...
$ teachers : num 10.9 11.1 82.9 14 71.5 ...
$ calworks : num 0.51 15.42 55.03 36.48 33.11 ...
$ lunch : num 2.04 47.92 76.32 77.05 78.43 ...
$ computer : num 67 101 169 85 171 25 28 66 35 0 ...
$ expenditure: num 6385 5099 5502 7102 5236 ...
$ income : num 22.69 9.82 8.98 8.98 9.08 ...
$ english : num 0 4.58 30 0 13.86 ...
$ read : num 692 660 636 652 642 ...
$ math : num 690 662 651 644 640 ...
# 下記いくつの データサイエンス パッケ一ジを library()関数で読み出します
# data science packages
library(tidyverse)
library(dplyr)
library(magrittr)# 少人数教育を測る変数(st.ratio: student-teacher ratio)を作る
# compute STR and append it to CASchools
# dplyr::mutate
CASchools %<>% mutate(st.ratio = students/teachers)# クラスサイズ変数の構築: d
# construct the class-size variable: d
# d =small_class, STR <20; d =large_class, STR>=20
CASchools %<>% mutate(d =rep("large_class", length(CASchools$district)))
CASchools$d[CASchools$st.ratio < 20] <- "small_class"
# dをカテゴリ変数にする
# make variable d into a categorical variable
CASchools$d %<>% as.factor# 学生のパフォーマンスを測定するアウトカム変数 `score` を作成する
# score = 読解力と数学のテストの平均点です
# compute TestScore and append it to CASchools
CASchools %<>% mutate(score = (read+math)/2)# dataframe structure
str(CASchools)'data.frame': 420 obs. of 17 variables:
$ district : chr "75119" "61499" "61549" "61457" ...
$ school : chr "Sunol Glen Unified" "Manzanita Elementary" "Thermalito Union Elementary" "Golden Feather Union Elementary" ...
$ county : Factor w/ 45 levels "Alameda","Butte",..: 1 2 2 2 2 6 29 11 6 25 ...
$ grades : Factor w/ 2 levels "KK-06","KK-08": 2 2 2 2 2 2 2 2 2 1 ...
$ students : num 195 240 1550 243 1335 ...
$ teachers : num 10.9 11.1 82.9 14 71.5 ...
$ calworks : num 0.51 15.42 55.03 36.48 33.11 ...
$ lunch : num 2.04 47.92 76.32 77.05 78.43 ...
$ computer : num 67 101 169 85 171 25 28 66 35 0 ...
$ expenditure: num 6385 5099 5502 7102 5236 ...
$ income : num 22.69 9.82 8.98 8.98 9.08 ...
$ english : num 0 4.58 30 0 13.86 ...
$ read : num 692 660 636 652 642 ...
$ math : num 690 662 651 644 640 ...
$ st.ratio : num 17.9 21.5 18.7 17.4 18.7 ...
$ d : Factor w/ 2 levels "large_class",..: 2 1 2 2 2 1 2 1 2 1 ...
$ score : num 691 661 644 648 641 ...
# summary statistics
summary(CASchools) district school county grades
Length:420 Length:420 Sonoma : 29 KK-06: 61
Class :character Class :character Kern : 27 KK-08:359
Mode :character Mode :character Los Angeles: 27
Tulare : 24
San Diego : 21
Santa Clara: 20
(Other) :272
students teachers calworks lunch
Min. : 81.0 Min. : 4.85 Min. : 0.000 Min. : 0.00
1st Qu.: 379.0 1st Qu.: 19.66 1st Qu.: 4.395 1st Qu.: 23.28
Median : 950.5 Median : 48.56 Median :10.520 Median : 41.75
Mean : 2628.8 Mean : 129.07 Mean :13.246 Mean : 44.71
3rd Qu.: 3008.0 3rd Qu.: 146.35 3rd Qu.:18.981 3rd Qu.: 66.86
Max. :27176.0 Max. :1429.00 Max. :78.994 Max. :100.00
computer expenditure income english
Min. : 0.0 Min. :3926 Min. : 5.335 Min. : 0.000
1st Qu.: 46.0 1st Qu.:4906 1st Qu.:10.639 1st Qu.: 1.941
Median : 117.5 Median :5215 Median :13.728 Median : 8.778
Mean : 303.4 Mean :5312 Mean :15.317 Mean :15.768
3rd Qu.: 375.2 3rd Qu.:5601 3rd Qu.:17.629 3rd Qu.:22.970
Max. :3324.0 Max. :7712 Max. :55.328 Max. :85.540
read math st.ratio d
Min. :604.5 Min. :605.4 Min. :14.00 large_class:181
1st Qu.:640.4 1st Qu.:639.4 1st Qu.:18.58 small_class:239
Median :655.8 Median :652.5 Median :19.72
Mean :655.0 Mean :653.3 Mean :19.64
3rd Qu.:668.7 3rd Qu.:665.9 3rd Qu.:20.87
Max. :704.0 Max. :709.5 Max. :25.80
score
Min. :605.5
1st Qu.:640.0
Median :654.5
Mean :654.2
3rd Qu.:666.7
Max. :706.8
# 数値変数を引き出し、その相関係数を計算します。
# これは、どのような変数を回帰式に含めるべきかのヒントを与えてくれる。
# 次のようなデータ可視化により、これらの変数間の関係を簡単に識別することができる。
# (Pearson) linear correlation: outcome, target, controls
target_controls <- select(CASchools, score, st.ratio, calworks:english)
cor_matrix <- cor(target_controls)
# Kendall coefficient: dependence
# cor(target_controls, method = c("kendall"))
# correlation matrix plot
# install.packages("corrplot")
library(corrplot)
corrplot::corrplot(cor_matrix, method = "circle")#corrplot::corrplot(cor_matrix, method = "color", addCoef.col = "grey")研究課題:少人数教育 (st.ratio) が学生の成績 (score) に与える定量的効果は何か?
# キレイなグラフを作るためのパッケージ
# data visualization, # qqplot2::qqplot
library(ggplot2)# scatter plot
#install.packages(c("ggExtra", "ggthemes"))
library(ggExtra)
library(ggthemes)
p <- ggplot(data =CASchools, aes(y =score, x =st.ratio)) +
geom_point(alpha = .5) +
theme_tufte(ticks=F)
ggMarginal(p, type = "histogram", fill="transparent")データ可視化を上手に行うためのお役立ちリンク集 (Useful links for doing good data visualization)
Rによる回帰分析 (Regression in R)
単回帰分析 (Simple regression illustrated through interactive graphs)
回帰式
\[ \hat{y}_i = \hat{\beta}_0 + \hat{\beta}_1 x_i \]
lm(y ~ x)
回帰分析は lm() を利用することで実行できます。
Rで回帰分析を行うと、推定されたバラメ一タ \((\hat{\beta}_0, \ \hat{\beta}_1)\)
の標準誤差とともに、有意差検定 (t検定) の結果を自動で出力してくれます。
回帰分析の結果を保存したオブジェクトを summary()
で要約すると、標準誤差は standard.error
としてレポ一トされており、有意差検定の結果は
t-statistics、p-value
としてそれぞれの変数に対して表示されています。
処置変数 d (カテゴリ変数)
\(d=1\) if
st.ratio<20; otherwise, \(d=0\).
回帰式
\[ \hat{score}_i = \hat{\beta}_0 + \hat{\beta}_1 d_i \]
# regress score on (d: discrete r.v.)
CASchools %$% lm(score ~ (st.ratio<20)) %>% summary
Call:
lm(formula = score ~ (st.ratio < 20))
Residuals:
Min 1Q Median 3Q Max
-50.496 -14.029 -0.346 12.884 49.504
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 650.077 1.393 466.666 < 2e-16 ***
st.ratio < 20TRUE 7.169 1.847 3.882 0.00012 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 18.74 on 418 degrees of freedom
Multiple R-squared: 0.0348, Adjusted R-squared: 0.0325
F-statistic: 15.07 on 1 and 418 DF, p-value: 0.0001202
クイズ: \(\hat{\beta}_0\) と \(\hat{\beta}_1\) の解釈は?ゼロとは有意に違うのでしょうか?
# load the color palettes
library(wesanderson)
# names(wes_palettes)
# e.g., Moonrise2, GrandBudapest1, GrandBudapest2# グラフ
# small class size (st.ratio < 20) v.s. large calss size
# scatter plots, linear specification
sp <- CASchools %>%
ggplot(aes(y =score, x =st.ratio, color =d)) +
geom_point(alpha = .7) +
stat_smooth(method = "lm", formula = y ~ x)
# d is a facotr variable with 2 levels; therefore n=2
sp + scale_color_manual(values=wes_palette(n=2, name="Moonrise2"))処置変数 st.ratio (連続変数)
回帰式
\[ \hat{score}_i = \hat{\beta}_0 + \hat{\beta}_1 st.ratio_i \]
# regress score on (st.ratio: continuous r.v.)
lm.fit.1 <- CASchools %$% lm(score ~ st.ratio)
summary(lm.fit.1)
Call:
lm(formula = score ~ st.ratio)
Residuals:
Min 1Q Median 3Q Max
-47.727 -14.251 0.483 12.822 48.540
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 698.9329 9.4675 73.825 < 2e-16 ***
st.ratio -2.2798 0.4798 -4.751 2.78e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 18.58 on 418 degrees of freedom
Multiple R-squared: 0.05124, Adjusted R-squared: 0.04897
F-statistic: 22.58 on 1 and 418 DF, p-value: 2.783e-06
クイズ: \(\hat{\beta}_0\) と \(\hat{\beta}_1\) の解釈は?ゼロとは有意に違うのでしょうか?
# グラフ
# 網掛け部分は95%信頼区間を表します
# scatter plot, linear specification
CASchools %>%
ggplot(aes(y =score, x =st.ratio)) +
geom_point(alpha = .5) +
stat_smooth(method = "lm", formula = y ~ x, color='#D55E00')重回帰分析 (Multiple regression)
次の2つのグラフは、なぜ共変量 (control variable) english
を含めるべきかを示しています。
したがって、重回帰分析を利用することで、より選択バイアスの影響が少ない分析ができます。
# motivation: multiple regressions, 2D scatter plot
# english: the share of English learning students
CASchools %>%
ggplot(aes(y =score, x =st.ratio, color =english)) +
geom_point(alpha = .5) # motivation: multiple regressions, 3D scatter plot
# english: the share of English learning students
#install.packages("plotly")
library(plotly)plot_interactive <- plot_ly(CASchools, x = ~st.ratio, y = ~english, z = ~score,
marker = list(color = ~english, colorscale = c('#FFE1A1', '#683531'), showscale = TRUE)) %>%
add_markers() %>%
layout(scene = list(xaxis = list(title = 'st.ratio'),
yaxis = list(title = 'PctEnglish'),
zaxis = list(title = 'score')),
annotations = list(
x = 1.13,
y = 1.05,
text = 'share of English learning students',
xref = 'paper',
yref = 'paper',
showarrow = FALSE
))
plot_interactive回帰式
\[ \hat{y}_i = \hat{\beta}_0 + \hat{\beta}_1 x_{1i} + \hat{\beta}_2 x_{2i} \]
# 重回帰分析
# multiple regression
lm.fit.2 <- CASchools %$% lm(score ~ st.ratio + english)
summary(lm.fit.2)
Call:
lm(formula = score ~ st.ratio + english)
Residuals:
Min 1Q Median 3Q Max
-48.845 -10.240 -0.308 9.815 43.461
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 686.03224 7.41131 92.566 < 2e-16 ***
st.ratio -1.10130 0.38028 -2.896 0.00398 **
english -0.64978 0.03934 -16.516 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 14.46 on 417 degrees of freedom
Multiple R-squared: 0.4264, Adjusted R-squared: 0.4237
F-statistic: 155 on 2 and 417 DF, p-value: < 2.2e-16
クイズ: なぜ、
st.ratioの推定効果が -2.2 (単回帰) から -1.1 (重回帰) に減少します?
不均一分散ロバスト標準誤差 (Heteroskedasticity-robust standard errors)
不均一分散とは?
テストの得点の分散が、生徒と教師の比率とともに増加することが観察されます。
estimatrは、一般的によく使われる線形推定量を提供するRパッケージであり、速度と使いやすさを考慮して設計されています。ユーザは、クラスタランダム化、ブロックランダム化、ブロックおよびクラスタランダム化を反映した線形推定量を選択することができる。
estimatr is an R package providing a range of
commonly-used linear estimators, designed for speed and for ease-of-use.
Users can choose an estimator to reflect cluster-randomized,
block-randomized, and block-and-cluster-randomized designs. https://github.com/DeclareDesign/estimatr
se_type: (ロバスト標準誤差の設定)
- classical: homoskedasticity (均一分散標準誤差,
lm()関数で使用される既定の設定です) - HC0: \(\omega_i = \hat{u}_i^2\). White (1980)
- stata: same as
vec(robust)used by Stata - HC1: \(\omega_i = \frac{n}{n-k}\hat{u}_i^2\)
- HC2: \(\omega_i=\frac{1}{1-h_i}\hat{u}_i^2\)
- HC3: \(\omega_i=\frac{1}{(1-h_i)^2}\hat{u}_i^2\)
lm_robust(y ~ x, se_type = "HC2")
library(estimatr)
# estimatr::lm_robust
lm.fit.2.robust <- CASchools %$%
lm_robust(score ~ st.ratio + english, se_type = "HC2")
summary(lm.fit.2.robust)
Call:
lm_robust(formula = score ~ st.ratio + english, se_type = "HC2")
Standard error type: HC2
Coefficients:
Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF
(Intercept) 686.0322 8.75433 78.365 1.278e-251 668.8241 703.2404 417
st.ratio -1.1013 0.43417 -2.537 1.156e-02 -1.9547 -0.2479 417
english -0.6498 0.03111 -20.888 7.919e-67 -0.7109 -0.5886 417
Multiple R-squared: 0.4264 , Adjusted R-squared: 0.4237
F-statistic: 222.8 on 2 and 417 DF, p-value: < 2.2e-16
st.ratioに対応する不均一分散ロバスト標準誤差は 0.43417
であり、これは均一分散の下で計算された標準誤差 (0.38028)
より大きいです。
ノンパラメトリック・ブートストラップ (Nonparametric bootstrap)
ブートストラップは、サンプリング分布を構築するための計算アルゴリズムです。The bootstrap is a computational algorithm for constructing sampling distributions.
- 理論的な近似に依存するのではなく、現在のサンプル(経験的データ分布)から再サンプリングを行うことで、サンプリング分布を模倣します。Instead
of relying on theoretical approximation, the bootstrap uses resampling
from your current sample (
empirical data distribution) to mimic the sampling distribution.
ALGORITHM 1. The Nonparametric Boostrap
Given data \(\{z_i \}_{i=1}^n := \{y_i, x_i \}_{i=1}^n\), for \(b =1, \ldots, B\):
Resample
with-replacement\(n\) observations \(\{z_i^b \}_{i=1}^n\)Calculate your estimate, say, \(\hat{\beta}_b\), using this resampled set.
Then \(\left\{\hat{\beta}_b \right\}_{b=1}^B\) is an approximation to the sampling distribution for \(\hat{\beta}\).
The standard deviation of your sampling distribution can be approximated by the standard deviation of this bootstrap sample: \[ \text{se}(\hat{\beta}) \approx \text{sd}(\hat{\beta}_b) = \sqrt{\frac{1}{B}\sum_b\left(\hat{\beta}_b -\hat{\beta}\right)^2}. \]
So long as the estimator is
unbiased, which means that \(\mathbb{E}[\hat{\beta}]=\beta\), then we can use this standard error to build the usual \(95\%\) confidence interval: \[ \beta\in \hat{\beta} \pm 2\text{sd}(\hat{\beta}_b). \]
The folloiwng algorithm applies to the case where \(\mathbb{E}[\hat{\beta}]\neq\beta\).
ALGORITHM 2. Nonparametric Bootstrap for Confidence Intervals
Given data \(\{z_i \}_{i=1}^n\) and full sample estimate \(\hat{\beta}\) for parameter \(\beta\), for \(b=1, \ldots, B\):
Resample
with-replacement\(n\) observation \(\{z_i^b \}_{i=1}^n\).Calculate your estimate, say, \(\hat{\beta}_b\), using this resampled set.
Calculate the error, say, \(e_b = \hat{\beta}_b - \hat{\beta}\).
Then \(\{e_b\}_{b=1}^B\) is an approximation to the distribution of errors between estimates and their targe.
- To get the \(95\%\) confidence interval for the true \(\beta\), we calculate the \(2.5\)th and \(97.5\)th percentiles of \(\{e_b\}_{b=1}^B\) – say, \(t_{0.025}\) and \(t_{0.975}\) – and set our interval as \[ \left[\hat{\beta} -t_{0.975}, \ \hat{\beta} -t_{0.025} \right]. \]
ブートストラップが破綻するのはいつか? When does the bootstrap break?
ブートストラップが破綻するのは、経験的データ分布が母集団の分布を適切に近似できていない場合や、対象とする統計量が母集団全体の分布に関する多くの情報を必要とする場合です。If the empirical data distribution is a poor approximation to the population or if the statistic you are targeting requires a lot of information about the entire distribution.
例えば、高次元のパラメーターの場合、母集団内の全ての変数間の依存関係を含む結合サンプリング分布を考慮する必要があるため、ブートストラップを信用しない傾向があります。e.g., we tend to not trust bootstrap for high-dimensional parameters because the joint sampling distribution is a function of all the population cross-variable dependencies.
このような場合には、ブートストラップの代替手法を検討する必要があります。When this happens, we look to alternative versions of the bootstrap.
Bootstrap in action, Take 1
B <- 1000
betas <- c()
for (b in 1:B) {
samp_b <- sample.int(nrow(CASchools), replace=TRUE)
reg_b <- lm(score ~ st.ratio + english, data=CASchools[samp_b,])
betas <- rbind(betas, coef(reg_b))
}
# note that betas is an array
# display the first five \hat{\beta}_b
betas[1:5,] (Intercept) st.ratio english
[1,] 685.2681 -1.083098 -0.6348671
[2,] 694.6744 -1.476707 -0.6560710
[3,] 687.2626 -1.133849 -0.6610360
[4,] 686.2406 -1.166699 -0.6136197
[5,] 685.7927 -1.108954 -0.6408945
# calculate the standard deviation of this boostrap sample
# apply(matrix, dimcode, function)
apply(betas[, 1:3], 2, sd)(Intercept) st.ratio english
8.88565009 0.44120502 0.03180366
Bootstrap in action, Take 2
# Nonparametric bootstrapping
library(boot)
# create a function producing beta_hat based on resampled observations
# this user-defined function will be used by boot::boot
beta_hat <- function(formula, data, indices) {
d.resampled <- data[indices,] # allows boot to select sample
fit <- lm(formula, data = d.resampled)
return(coef(fit))
}
# regression specification
outcome <- "score"
covariates <- c("st.ratio", "lunch", "income", "english")
reg_spec <- as.formula(
paste(outcome,
paste(covariates, collapse = " + "),
sep = " ~ "))
print(reg_spec)score ~ st.ratio + lunch + income + english
# bootstrapping with 1000 replications
boot_results <- boot(data = CASchools, statistic = beta_hat,
formula = reg_spec, R=1000)
# standard errors - bootstrap, 1 intercept + 4 regressors
# index=1: intercept; index=2: st.ratio; ...
boot_results
ORDINARY NONPARAMETRIC BOOTSTRAP
Call:
boot(data = CASchools, statistic = beta_hat, R = 1000, formula = reg_spec)
Bootstrap Statistics :
original bias std. error
t1* 675.6082083 -1.083597e-01 6.03556103
t2* -0.5603890 2.180953e-03 0.24953695
t3* -0.3963660 1.280183e-04 0.03056716
t4* 0.6749840 2.448742e-03 0.08257985
t5* -0.1943284 5.797645e-05 0.03397513
# the bootstrap sample mean of beta_hat_st.ratio = -0.5582080
colMeans(boot_results$t)[1] 675.4998487 -0.5582080 -0.3962380 0.6774327 -0.1942704
# thus the bias = the bootstrap sample mean - the estimate based the original sample
# -0.5582080 - (-0.5603890) = 0.002180953
# when constructing CI, remember to do bias-correction# histogram of beta_hat corresponding to the st.ratio variable
# index=1: intercept; index=2: st.ratio; ...
plot(boot_results, index=2)# get 95% confidence interval of beta_hat corresponding to st.ratio
boot.ci(boot_results, conf = c(0.90, 0.95), type = "perc", index=2) #st.ratioBOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 1000 bootstrap replicates
CALL :
boot.ci(boot.out = boot_results, conf = c(0.9, 0.95), type = "perc",
index = 2)
Intervals :
Level Percentile
90% (-0.9808, -0.1819 )
95% (-1.0714, -0.0903 )
Calculations and Intervals on Original Scale
偽発見率の制御 (False Discovery Rate Control)
多重性の問題。回帰推定の問題を考えます。100個の係数のうち、5つだけが結果と実際に関連しているとします。The problem of multiplicity. Imagine a regression estimation problem where \(5\) out of \(100\) coefficients have a real relationship with the result.
この場合、真の係数セットはスパース(疎)であり、多くの\(\beta_j\)がゼロであると言えます。We would say here that the true coefficient set is sparse: most of the \(\beta_j\) are zero.
最良のケース では、仮説検定を実行して、これら5つの真のシグナルをすべて発見したと仮定します(すなわち、偽陰性はない)。In the best-case scenario, when you run hypothesis tests, you find all of these five true signals (i.e., there are no false negatives).
ここで、残りの\(95\)個の係数に対して、\(\alpha=0.05\)の有意水準を用いて検定を行うとします。この場合、有用でない95個の変数の\(5\%\)について誤って有意であると結論することが期待されます。これにより、最終モデルに4.75(約5)の偽の説明変数(スプリアス・リグレッサー)が含まれることになります。Let’s assume you do find all five and test the remaining \(95\) with \(\alpha=0.05\) significance cutoff. In that case, you expect that you will erroneously conclude significance for \(5\%\) of the useless \(95\) variables, leading you to include \(4.75\approx 5\) spurious regressors in the final model.
したがって、False Discovery Proportion (FDP)(偽発見率)は次のように定義されます:Thus the false discovery propotion (FDP) \[ FDP = \frac{\# \ \text{false positives}}{\# \ \text{tests called significant}} = \frac{5}{10} = 50\% \]
FDP(False Discovery Proportion)は、適合したモデルの性質です。これを正確に知ることはできません。つまり、誤りを犯しているかどうかは不明です。しかし、その期待値である偽発見率 (FDR: False Discovery Rate) を制御することは可能です。The FDP is a property of the fitted model. You cannot know it: either you have made mistakes or you have not. However, you can control its expectation, the false discovery rate (FDR): \[ FDR = \mathbb{E}[FDP] =\mathbb{E}\left[\frac{\# \ \text{false positives}}{\# \ \text{tests called significant}} \right]. \]
次のアルゴリズムを用いることで、選択した \(q\)-カットオフ(例:一般的には \(0.1\))に対して \(FDR \leq q\) を確実に保証することができます。You can ensure that \(FDR \leq q\) for some chosen \(q\)-cutoff (e.g., \(0.1\) is commom) in light of the following algorithm.
ALGORITHM 3. Benjamini-Hochberg (BH) FDR Control
For \(N\) tests, with \(p\)-values \(\{p_1, \cdots, p_N \}\) and target FDR \(q\):
Order your \(p\)-values from smallest to largest as \(p_{(1)}, \cdots, p_{(N)}\).
Set the \(p\)-value cutoff as \(p^*=\max\left\{p_{(k)}: p_{(k)}\leq q\frac{k}{N}\right\}\).
The rejection region is then the set of all \(p\)-values \(\leq p^*\), and this ensure that \(FDR\leq q\). \(\quad \square\)
このアルゴリズムは、検定間の独立性を仮定していることに注意してください。Notice that the algorithm assumes independence between tests.
FDR control in action
# multiple regressions: 7 covariates
lm.fit.7 <- CASchools %$% lm_robust(score ~ st.ratio + calworks + lunch + computer +
expenditure + income + english, se_type = "HC2")
summary(lm.fit.7)$coef Estimate Std. Error t value Pr(>|t|) CI Lower
(Intercept) 6.618848e+02 1.011465e+01 65.4382631 8.746942e-220 6.420021e+02
st.ratio -2.782033e-01 3.131453e-01 -0.8884160 3.748354e-01 -8.937650e-01
calworks -9.174277e-02 6.369565e-02 -1.4403303 1.505334e-01 -2.169518e-01
lunch -3.709982e-01 4.470528e-02 -8.2987567 1.518469e-15 -4.588771e-01
computer 4.100852e-04 8.207297e-04 0.4996593 6.175817e-01 -1.203255e-03
expenditure 1.743849e-03 9.054684e-04 1.9259079 5.480367e-02 -3.606548e-05
income 6.190124e-01 8.463125e-02 7.3142304 1.361586e-12 4.526495e-01
english -2.113803e-01 3.889648e-02 -5.4344319 9.427272e-08 -2.878406e-01
CI Upper DF
(Intercept) 681.767588421 412
st.ratio 0.337358497 412
calworks 0.033466221 412
lunch -0.283119321 412
computer 0.002023425 412
expenditure 0.003523763 412
income 0.785375360 412
english -0.134919948 412
# Benjamini-Hochberg (BH) False Discovery Rate Control,
# Taddy (2019), p. 32, Algorithm 3 ensuring FDR \leq q, say, 0.1
# extract the p-values for the nine regression coefficients (excluding the intercept)
# and then plot p-values against their rank from smallest to largest
# slope of the BH cutoff line: q/N, where N is sample size
# decision rule: find the largest p-value that is below this line
# and call it and all smaller p-values "significant,
# so that you expect that aroud q = 0.1 (FDR) of them are false positives
# grab p-values
#-c(1,2): drop the intercept and st.ratio
# since st.ratio is the target variable, it should remain in the model specification
pvals <- summary(lm.fit.7)$coef[-c(1,2), "Pr(>|t|)"]
# plot them: it looks like we have some signal here
hist(pvals, xlab="p-value", main="", col="lightblue")# cutoff function
BH_cutoff <- function(pvals, q=0.1){
pvals <- sort(pvals[!is.na(pvals)])
N <- length(pvals)
k <- rank(pvals, ties.method="min")
alpha <- max(pvals[ pvals<= (q*k/N) ])
plot(pvals, log="xy", xlab="order", main=sprintf("FDR of %g",q),
ylab="p-value", bty="n", col=c(8,2)[(pvals<=alpha) + 1], pch=20)
#log="xy": both axes are to be logarithmic
#bty="n"": no box
lines(1:N, q*(1:N)/(N+1))
return(alpha)
}# Benjamini-Hochberg (BH) False Discovery Rate Control cutoff value
BH_cutoff(pvals) # 0.05480367[1] 0.05480367
# at 10% FDR, we get 3 `signif'
signif <- which(pvals <= 0.05480367) # multiple regressions: 3 covariates
ols.fit.4 <- lm_robust(score ~.,
data=CASchools[,c("score", "st.ratio", names(signif))], se_type = "HC2")
summary(ols.fit.4)
Call:
lm_robust(formula = score ~ ., data = CASchools[, c("score",
"st.ratio", names(signif))], se_type = "HC2")
Standard error type: HC2
Coefficients:
Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF
(Intercept) 675.6082 6.22466 108.537 9.114e-307 663.3724 687.84400 415
st.ratio -0.5604 0.25603 -2.189 2.917e-02 -1.0637 -0.05711 415
lunch -0.3964 0.03033 -13.066 6.385e-33 -0.4560 -0.33674 415
income 0.6750 0.08448 7.989 1.352e-14 0.5089 0.84106 415
english -0.1943 0.03336 -5.825 1.147e-08 -0.2599 -0.12875 415
Multiple R-squared: 0.8053 , Adjusted R-squared: 0.8034
F-statistic: 463.2 on 4 and 415 DF, p-value: < 2.2e-16
不良コントロール (Bad Controls)
回帰に職業ダミー(occupation dummy)をコントロール変数として含めるべきか?
以下の説明は、[Paul Hünermund氏のブログ記事]から抜粋したものです。
The following explaination is extracted from Paul Hünermund’s blogpost.
良いコントロ一ル \(X\) は、処置変数とアウトカム変数の両方に影響を与えますが、それをコントロールすれば(それを測定することができれば)問題ありません。観察されていない交絡因子 \(U\) の存在は、この場合問題ではない。
Good Controls. \(X\) affects both the treatment and outcome variable, but once we control for it (given that we can measure it) we’re fine – unconfoundedness holds. The presence of an unobserved confounder \(U\) doesn’t matter in this case.
不良コントロ一ル \(T\)と\(X\)の間の因果関係が、今では逆の方向に向かっているということです。大学教育 \((T)\) は、将来の職業 \((X)\) に影響を与え、その結果、将来の収入 \((Y)\) に影響を与えます。もし\(X\)をコントロールすると、問題が発生します。それは、\(U\)を通る因果関係の道を広げてしまう(\(T \rightarrow X \leftarrow U \rightarrow Y\), なぜなら、\(X\)はその道のコライダーだから)、バイアスがかかってしまうのです。
Bad Controls. Only a tiny detail has changed, namely that the causal link between \(T\) and \(X\) now goes in the other direction. College education \((T)\) affects your future occupation \((X)\), which in turn affects your future income \((Y)\). Suddenly, controlling for \(X\) causes problems. It opens up the causal path going through \(U\) (\(T \rightarrow X \leftarrow U \rightarrow Y\), because \(X\) is a collider on that path) and produces bias.
リーディング (Readings)
Angrist and Pischke (2015). Mastering ’Metrics: The Path from Cause to Effect. Princeton University Press.
西山慶彦、新谷元嗣、川口大司、奥井亮 (2019),『計量経済学』、有斐閣:
第4章 - 線形単回帰モデルの推定と検定 (pages 103 – 139)
第5章 - 重回帰モデルの推定と検定 (pages 141 – 207)
実証練習問題 (Empirical Exercises)
[E1]
GoogleClassroom には、デ一タファイル TeachingRatings
があります。 このデ一タはテキサス大学オ一スティン校の Daniel Hamermesh
教授から提供されたもので、Amy Parker 氏との共同研究 “Beauty in the
Classroom: Instructor’ Pulchritude and Putative Pedagogical
Productivity,” Economics of Education Review, August 2005,
24(4): pp. 369-376 で使用されたものです。
そこには、テキサス大学オ一スティン校での 463
の講義に関する授業評価、授業の特徴、教授の特徴が含まれています。
デ一タの詳細は、TeachingRatings_Description
に書かれています。 このデ一タセットの特徴の一つは、教授の「(外見上の)
美しさ」が 6
人の学生からなるパネルによって評価され、指標化されている点です。
この練習問題では、授業評価が教授の「美しさ」とどのように関係しているかを調べます。
The data file TeachingRatings contains data on course
evaluations, course characteristics, and professor characteristics for
463 courses at the University of Texas at Austin. A detailed description
is given in TeachingRatings_Description. One of the
characteristics is an index of the professor’s “beauty” as rated by a
panel of six judgers. In this exercise you will investigate how course
evaluations are related to the professor’s beauty.
教授の美しさ
Beautyと 平均授業評価Course_Evalとの散布図を作成しなさい。 2つの変間に関係があるように見えるでしょうか。 – Construct a scatterplot of average course evaluationsCourse_Evalon the professor’s beautyBeauty. Does there appear to be a relationship between variables?平均授業評価
Course_Evalを教授の美しさBeautyで回帰しなさい。 定数項の推定値、傾きの推定値はそれぞれいくつですか。 なぜ定数項の推定値はCourse_Evalの標本平均と等しいのでしょうか。 [ヒント: Beauty の標本平均はいくですか。] – Run a regression of average course evaluationsCourse_Evalon the professor’s beautyBeauty. What is the estimated intercept? What is the estimated slope? Explain why the estimated intercept is equal to the sample mean ofCourse_Eval.陳教授の Beauty の値は平均値であり、森教授の Beauty の値は、平均値より 1 標準偏差高い値です。 森教授と陳教授の授業評価を予測しなさい。 – Professor Chen has an average value of Beauty, while Proessor Mori’s value of Beauty is one standard deviation above the average. Predict Professor Mori’s and Professor Chen’s course evaluations.
回帰係数の大きさに関してコメントしなさい、
BeautyとCourse_Evalに及ぼす推定された効果は大きいですか、小さいですか、あなたが「大きい」あるいは「小さい」と評価する基準は何ですか。 この推定された効果が因果効果なのかどうかを説明してください。 – Comment on the size of the regression’s slope. Is the estimated effect ofBeautyandCourse_Evallarge or small? Explain what you mean by “large” and “small”. Is the estimated effect a causal effect? Explain.Beautyは授業間の評価の分散の大部分を説明していますか。 – DoesBeautyexplain a large fraction of the variance in evaluations across courses? Explain.授業の種類や教授の特徴をコントロ一ルするため、変数を追加して、
Course_EvalをBeautyで回帰しなさい。 具体的には、コントロ一ル変数にIntro、OneCredit、Female、Minority、NNEnglishを追加しなさい。BeautyがCourse_Evalへ及ぼす効果の推定値はいくつですか。 (2)の単回帰には、除外された変数のバイアスが含まれているでしょうか。 – In order to control for course characteristics and professor characteristics, we regressCourse_EvalonBeautywith several control variables in the regression specification includingIntro、OneCredit、Female、MinorityandNNEnglish. Based on this multiple regression, what is the effect ofBeautyonCourse_Eval? Does the simple regression in part (2) include the omitted variable bias?スミス教授は平均的な容姿の黒人男性で、英語を母国語としています。 彼は 3 単位もの上級科目を教えています。 スミス教授の授業評価を予測しなさい。 – Professor Smith is an average looking black male whose native language is English. He teaches 3 credits of advanced courses. Predict Professor Smith’s class evaluations.
[E2]
ここでは、アメリカにおける男女間の賃金格差を題材としたケーススタディを検討します。
Codebook.txt には労働者の職業関連特性の説明が記載されています。
Pay.discrimination.Rdataは、賃金および経験(exp)、性別、教育などの職業関連特性に関するCPS(2012年)のデータです。Here
we consider the case study of the gender wage gap in the United States.
Codebook.txt contains the description of worker job-relevant
characteristics. Pay.discrimination.Rdata: the CPS (2012)
data on wages and job-relevant worker characteristics, such as
experience (exp), gender, education.
同じ特性を持つ男性と女性の間で、賃金にどのような違いがあるでしょうか?回帰モデルの仕様をどのように決定したか説明してください。What is the difference in wages between men and women with the same characteristics? Explain how you determined the regression specification.
回帰推定量の異質分散性に頑健な標準誤差を推定してください。また、各回帰係数に対する仮説検定を実施してください。Estimate the heteroskedasticity-robust standard error of your regression estimator. Conduct hypothesis testing on each regression coefficient.
ブートストラップを使用して、回帰推定量の標準誤差を推定してください。Estimate the standard error of your regression estimator by using bootstrapping.